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PACS 64 . 70 . kp - Ionic crystals 
PACS 68 . 65 . Ac - Multilayers 

PACS 73 . 20 . -r - Electron states at surfaces and interfaces 

Abstract -We study the ground state structure of electronic-like bilayers, where different phases 
compete upon changing the inter-layer separation or particle density. New series representa- 
tions with exceptional convergence properties are derived for the exact Coulombic energies under 
scrutiny. The complete phase transition scenario -including critical phenomena- can subsequently 
be worked out in detail, thereby unifying a rather scattered or contradictory body of literature, 
hitherto plagued by the inaccuracies inherent to long range interaction potentials. 
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The prediction by Wigner that strongly correlated 
charge carriers in a uniform compensating background 
could crystallize [I] , was first realized experimentally with 
electrons at the surface of liquid Helium [2J, which form 
a two-dimensional structure. Since then, the study of low 
dimensional electronic systems has shown no abating and 
in particular, the bilayer geometry singles out. It appears 
significantly richer than its monolayer counterpart and has 
been investigated in different settings : GaAs quantum 
wells [3] or other semiconductors [3], quantum dots [5], 
graphene [BJ, boron nitride [TJ, laser-cooled trapped ion 
plasmas [5], dusty plasmas [3] and colloids jTQl|TT| . In 
light of these applications, it is essential to understand 
the ground state features of Coulombic bilayers, start- 
ing with the classical limit. This problem has received 
significant attention in the last 20 years, as such [12 1G 
or supplemented with finite temperature analysis [17H19) . 
In addition, ground state ordering impinges on strong- 
coupling expansions describing counterintuitive yet ubiq- 
uitous electrostatic phenomena, such as like-charge attrac- 
tion [2"01 - |2"5] or charge reversal J23jr25] ■ This body of work 
has revealed the main features of ground state structure, 
but there exist surprising discrepancies in the literature, 
especially on the respective domains of existence of the dif- 
ferent phases possible. The reason lies in the long range 
nature of Coulombic interactions, a common bane for such 
analysis. Our goal here is to resolve existing controversies, 
precisely locate all phases and discuss the critical behav- 
ior associated. All results reported are exact; they are 
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Fig. 1: Side view of the two parallel plates at a separation 
d. From their common uniform surface charge density —aq, 
we define the dimensionless distance as n = dy/a. Ions, which 
bear charge q, are shown as black or white disks for visual ease, 
but they are point-like in the present study. 



obtained from new series representations of lattice sums 
for Coulomb law. 

We consider an ensemble of identical classical point 
charges q, interacting through a 1/r pair potential, and 
confined between two symmetric parallel charged walls. 
These boundaries both bear a uniform surface charge of 
density —<rq, so that global electroneutrality holds. At fi- 
nite temperature T, the charges do populate the interior of 
the slab. For T = though, the charges evenly condense 
on the opposing walls, thereby forming a bilayer ground 
state [26j[22] the structure of which depends on a single 
dimensionless parameter n = d^fa, where d is the inter- 
plate distance (see Fig. [JJ. It is known that five structures 
can be realized upon increasing 77; they will be referred to 
following standard terminology, common to the classical 
[TJJ and quantum contexts [55] . To begin with, the limits 
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Fig. 2: Schematic representation of the different ground states 
encountered when the dimensionless distance rj increases. The 
open and filled symbols show the locations of ions on the oppo- 
site surfaces (see Fig. [TJ. The arrows are for lattice vectors ai 
and a2, from which we define the aspect ratio A — |<i2|/|aij: 
A = v3 and 1 for structures I and III respectively, while struc- 
ture II interpolates between I and III with 1 < A < v^3- For 
structure IV, the order parameter is the angle tp between ai 
and 02. We have p = 7r/2 for structure III whereas p = n/3 
for structure V, and in general, tv/3 < tp < n/2. Structures 
I, III and V are rigid, as opposed to the soft cases II and IV 
where the unit cell geometry depends on inter-plate separation, 
through A and tp, respectively. Note that the shift between the 
two opposite crystals materialized by open and filled symbols 
is (ai + a,2)/2 for structures I, II, III, IV but (ai + a2)/3 for 
structure V. 



of small and large rj are both straightforward. For 77 — > 0, 
a genuine two-dimensional one component plasma is pro- 
duced [20] . where the strong mutual repulsion between 
charges leads to a triangular Wigner crystal [30], the so- 
called structure I. Conversely, for 77 — > 00, the two lay- 
ers decouple and a hexagonal crystal forms on each plate 
(structure V). These two crystals adopt a staggered con- 
figuration, to minimize inter-layer repulsion. For interme- 
diate reduced distances 77, three other structures are met, 
see Fig. [5J a staggered rectangular lattice (structure II), 
a staggered square lattice (structure III) , and a staggered 
rhombic lattice (structure IV). Note that while one can 
evolve continuously through the sequence I — > II — > III — > 
IV, no continuous deformation allows to create structure 
V from one of the others. The transitions between phases 
will therefore be of different orders, with characteristics 
and critical exponent (in the continuous cases) that will 
be worked out explicitly below. A goal of our analysis is 
to precisely locate the transition points between phases: 
indeed, a dispersion of about 20% exists for the hitherto 
reported threshold rjiy between structures IV and V, see 
13, 14, 16, 18:. In addition, controversial results have been 
reported for the transition point rji between structures I 
and II: 771 ~ 0.006 from Ewald summation technique [13] . 
771 ~ 0.028 from Monte Carlo simulations |T51 , whereas 
lattice sum minimization of Yukawa systems in the un- 



screened limit hints at 771 = |15j , meaning that structure 
I could possibly only exist at precisely 77 = 0. 

We start by addressing the question whether rj\ = 
or 7^ 0, and to this end, we compute the energy E(A, 77) 
of structure II. For a given layer, the 2D lattice points 
are indexed by ja\ + ka-i where j and k are integers and 
the lattice vectors a\ = a(l,0), = a(0, A) are shown 
in bold in Fig. [3] The global electroneutrality requires 
that a 2 a A = 1. The aspect ratio A fulfills 1 < A < \/3 
with A = \/3 for structure I and A = 1 for the (square) 
structure III. The dielectric constant of the medium is set 
to unity for the sake of simplicity, and the total energy 
per ion E(A,vj) is written as the sum of intra- and inter- 
layers contributions. We first restrict ourselves to a disk 
of finite radius R around a given reference ion located at 
(0,0). Considering ion- ion and ion-plate interactions, we 
have the intra-layer energy 
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with the restriction j 2 + k 2 A 2 < (R/a) 2 . It is expedient 
and common procedure [191131] to use the gamma identity 
(n/z) 1 /' 2 = J °° t~ 1 / 2 exp(—zt) dt valid for z > 0, which 
provides us with a simple expression where the limit R — > 
00 can be readily taken: 
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Here, the second line is obtained from the substitution 
tA — > t and the condition a 2 a A = 1; the third line stems 
from considering separately the domains t £ [0, tt] and 
t G [it, 00} in the integral, substituting ir 2 /t — > t and sub- 
sequently using Poisson summation formula 
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The inter-layer energy contribution -Einter is amenable to a 
similar treatment [32) . and the last step of the procedure 
consists in introducing the function 
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for 77 > 0. 
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We finally end up with the series representation for the 
total energy E(A, 77) = -Ejntra + Winter: The function 
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z,, generalizes to two-layer problems the so-called Misra 
function [33J used extensively in single-layer lattice sum- 
mation [3"fl!3o1 - Our use of §5§ will be three-pronged: 
it allows to show analytically that 771 = 0, to calcu- 
late explicitly the singular behavior near critical points 
and is moreover particularly suited for numerical eval- 
uations. From an operational point of view, the scries 
([5]) is indeed endowed with remarkable properties: it is 
free of singular terms, and importantly, converges ex- 
tremely quickly. The error made upon truncating the 
series in the energy expression ([5]) at order j = k = AI 
behaves like exp(— cM 2 )/M, where c is a constant of or- 
der unity. We first document the convergence property 
on the single-layer case of structure I, for which the exact 
energy is E(y/S, 0)/(q 2 V2^) = -1.96051578931989165 . . ., 
which is directly the Madclung constant of the 2D hexag- 
onal Wigncr crystal. Cutting the series ([5]) at M, we ob- 
tain the exact value with a precision of 2,5, 10, 17 digits 
with notably small cutoffs M = 1, 2, 3, 4 respectively. This 
makes numerical calculations extremely fast and efficient 
on any workstation. A similar accuracy is met for all other 
structures and parameter values reported here, and all nu- 
merical results quoted below have been obtained with the 
cutoff M = 5. 

We now turn our attention to the threshold 771 which 
defines the stability window of structure I. For a given 
distance 77, we proceed by calculating the Taylor expansion 
of ([5]) in the small parameter e = \/3 — A, which yields 



f 2 ( v ) = 0.0408440789 . . . + 0(j] 2 ). As a consequence, 
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where and the function f 2 is also explicitly known [32] , 
To investigate the stability of structure I, it is suffi- 
cient to study the sign of /1, which is worked out from 
a Taylor expansion for small 77. The first two deriva- 
tives of this function fi vanish at 77 = and we have 
/i(t7) = -0. 5833059875... r/ 2 + O^ 4 ), hence an energy 
gain upon increasing e as compared to the e = case 
(structure I). This implies that 771 = 0: at finite but 
small distances 77, the optimal phase is not structure I. 
To obtain the optimal value of e selected and that we 
denote e*, we further Taylor expand 72(77) which yields 



V3 - A* = e* 
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7.14064... 77 2 + C(?7 4 ), (8) 



which entails that the energy change scales like 77 4 . For 
the thresholds 771 reported in Refs. Q31[TH], a relative ac- 
curacy of 10 -9 was therefore required to answer the finite 
or vanishing 771 question. The accuracy of our findings is 
□ 




Fig. 3: The difference between the dimensionless energies 
[E(A,7]) - S(v / 3,7?)]/(g 2 v / 2^) versus e = V3 - A, for 77 = 
5. 10 -3 . The analytical formula ^ with the Taylor expansions 
of /1 and /2 given in the text is shown by the continuous line. 
It is compared to the numerical evaluation of the series ([5| with 
a cutoff M = 5 (symbols). The optimal e value following from 
the prediction © is shown by the dashed vertical line. 

The above analysis shows that the evolution from struc- 
ture I to structure II is not a phase transition in the com- 
mon sense. The situation differs between structures II and 
III. To inspect the corresponding transition, we note that 
E(A,rj) enjoys the symmetry A — > A -1 , as is clear from 
Fig. [5] where a global rotation of 7r/2 does not affect the 
energy but interchanges lattice vectors a\ and a 2 . The 
value A = 1 characterizing structure III is therefore a 
self-dual point, and it will now be convenient to param- 
eterize the aspect ratio as A = cxp(e). All expressions 
will then be even in e. The expansion of E(e e ,r)) in small 
e-deviations yields 



£(e £ ,77) E(1, V ) 
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Fig. 4: The transition II — >• III: Test of the analytical asymp- 
totic relation (|10|) (dashed line) against numerical minimization 
of the energy Q (solid curve), in a log- log scale. The numerical 
data of Ref. [14] are shown by the circles. 



where g 2 and 34 are explicitly known |32j , in a series form 
that is very reminiscent of Eq. (|7|). The bilayer energy 
appears in a standard Ginzburg-Landau form [35], but it 
should be emphasized that at variance with mean- field ar- 
guments usually underlying such approaches, our expres- 
sion is exact. The critical point 7711 sought for is the root 
of g 2 (r)) = 0, which gives 7711 = 0.2627602682... It ap- 
pears here that the thresholds reported in earlier works 
were accurate: 0.27 [H], 0.262 [H], 0.27 [H] and 0.28 [IB]. 
Proceedings along similar lines as for the I — > II crossover, 
we Taylor expand (72(77) and <74 (77) to leading order around 
7711. The former behaves like (77 — r)u) while the latter is 
constant, a prototypical scenario for a continuous phase 
transition with critical index (3 = 1/2 [35]. Specifically, 
we get 



A* 
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2g 4 (?7) 



1/2 



1. 48031 (10) 



This expression applies for 77 < 7711, in the stability do- 
main of structure II, and is in excellent agreement with 
our numerical calculations from Eq. ([5]), see Fig. [4] 

The task remaining is to find the series representations 
for structures IV and V. We first address structure IV. 
Implementation of the procedure that led to the series ([5j 
becomes possible once the distance between a reference ion 



and an arbitrary ion located on the same layer at r(j, k) 
jai + ka 2 is expressed as 



\r(j,k)\ 2 
= a 2 [(J - 



a 2 (j 2 + k 



2jk cos ipj 

kf cos 2 ((^/2) + (j - k) 2 sin 2 (^/2)] (11) 



The latter "diagonalized" form in terms of indices, pro- 
vides the starting point to write the intra-layer Coulomb 
energy (summing l/y/\r{j, k)\ 2 ), and suggests to intro- 
duce new indices n and m: if j + k is even, we define 
n = (j + k)/2 and m = (j — k)/2. If j + k is odd, we intro- 
duce indices n = (j + k + 1)/2 and n = (j — k + l)/2. Like- 
wise for the inter-layer interactions, taking due account of 
the shift (a 1 + a 2 )/2 between opposite layers. Building on 
the gamma identity and Poisson summation formula, the 
series form for the energy E\y ensues [35]. This energy 
depends on the angle tp and of course on the distance 77. 
For our purposes, rather than the lengthy explicit form, it 
is sufficient to report the Landau-like expansion of Ejy in 
the vicinity of <p = it/2. A convenient expansion param- 
eter is e such that exp(e) = tan(y>/2), and the invariance 
tp — > 7r — ip makes E\\ an even function of e. In the small 
e region of interest associated to the vicinity of 7r/2 for tp, 
we obtain an expansion up to order e 4 of the same form 
as (0). This teaches us that structure III is unstable for 
t? > 7?iii = 0.6214809246 . . ., to be compared to the thresh- 
olds 0.61 [T3J, 0.622 QI], 0.62 [T5J, 0.59 [IE] while structure 
IV was not considered in [T^]. We furthermore again ob- 
tain a second order phase transition with critical index 1/2 
and explicit order parameter close to the transition point 
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in excellent agreement with our numerical data. 

The transition IV — > V is discontinuous, which made 
its characterization more elusive in previous publications. 
Our method, though, is easily adapted to the geometry 
of structure V. The series representation for E\ (77) should 
be compared to £rv (</?*, 77) evaluated for the optimal dis- 
tortion angle (p*(rf). Requiring that Ey(j)) = Eiv(tp*,rf), 
we obtain the last 77-thrcshold that was remaining to be 
specified: 7/iv = 0.73242 . . . Previous investigations gave 
0.75 [UJ, 0.732 [H], 0.87 [TB] and 0.70 [TSJ. For 77 > rj w , 
structure V is energetically favorable. As a by-product of 
our analysis, we report the large distance behavior of the 
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interplate pressure P = —2<rdEy/dd: 



REFERENCES 



P = — 6ir(aq) 2 exp 



-in 



V 



V23 1 / 4 

in agreement with [T4j but at variance with [T31I22] 





second order 


first order 

1 






I \ 

\ 


^Nv^m ^iv 







0.263 


0.621 0.732 




/ 


II i 77/ 


i IV V 





Fig. 5: Sequence of structures encountered as a function of 
reduced inter-plate separation 77. The values reported for the 
different stability thresholds are rounded to the third digit. 

To summarize, we have derived series representations 
for the different Coulomb lattice sums pertaining to the 
ground state of classical bilayer systems. The derivation, 
worked out explicitly for the five different structures that 
were known to compete at vanishing temperature, results 
from a series of transformation rooted in the general the- 
ory of Jacobi 9 functions 02- The resulting series pro- 
vide the thresholds delimiting the domains of validity of 
the different phases, that were prone to some fluctuations 
in previous works. Figure [5] provides an overview of our 
main findings. We could in particular show that the sim- 
ple hexagonal structure I can only exist in the limiting 
case of a vanishing interplate distance, and is preempted 
by a buckled phase for all 7] 7^ 0. This is the scenario 
first reported in Ref . [TS] , which differs from several other 
studies that assigned a finite stability window to phase I. 
Whereas the evolution I — > II is not a phase transition, we 
could show that the continuous transitions II — > III and 
III — > IV have critical index (3 = 1/2. In addition, our 
series representation is endowed with exceptional conver- 
gence properties, providing typically more than 10 digits 
of accuracy when retaining only the first 4 or 5 terms 
involved. Relinquishing the symmetry between the two 
plates to address the cases where they bear different sur- 
face charges is an interesting venue for future work. This 
brings the difficulty that local electroneutrality no longer 
holds at the single-plate level in the ground state [37] , ex- 
cept presumably at large separations. Our approach can 
also be extended to bilayers and multilayers with repulsive 
Yukawa or inverse-power-law interactions, that deserve at- 
tention. 
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